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ABSTRACT 

The propagation of a shock wave into an interstellar medium is investigated by two-dimensional nu¬ 
merical hydrodynamic calculation with cooling, heating and thermal conduction. We present results of 
the high-resolution two-dimensional calculations to follow the fragmentation due to the thermal insta¬ 
bility in a shock-compressed layer. We find that geometrically thin cooling layer behind the shock front 
fragments into small cloudlets. The cloudlets have supersonic velocity dispersion in the warm neutral 
medium in which the fragments are embedded as cold condensations. The fragments tend to coalesce 
and become larger clouds. 

Subject headings: ISM: clouds — ISM: molecules — ISM: supernova remnants — shock waves — 
turbulence 


1. INTRODUCTION 

The interstellar medium (ISM) and cold gas clouds are 
characterized by a clumpy substructure and a turbulent 
velocity field (Larson 1981). The maintenance and dis¬ 
sipation processes of the turbulence are supposed to be 
important in the theory of star formation (McKee 1989, 
Nakano 1998, Elmegreen 1999). The understanding of the 
origin of cold clouds and their internal substructure has 
therefore fundamental importance for a consistent theory 
of star formation and ISM. 

Some of the polarization maps and direct measurements 
of field strength in some star forming regions suggested the 
importance of MHD (Alfven waves) turbulence (Myers & 
Goodman 1988, Crutcher et al. 1993). Recent simula¬ 
tions of MHD turbulence, however, suggest that it dis¬ 
sipates rapidly (Gammie & Ostriker 1996, Mac Low et 
al. 1998, Mac Low 1999, Ostriker et al. 1999). Possi¬ 
ble sources of energy supply are winds and outflows from 
young stellar objects (Franco & Cox 1983, McKee 1989). 
Note however that clumpy structures with supersonic ve¬ 
locity dispersions are also observed even in regions where 
star formation is inactive. Thus, the origin of clumpy cloud 
structure remains as an outstanding issue. 

We propose that the clumpiness in clouds arises nat¬ 
urally from their formation through thermal instability 
which acts on timescales that can be much shorter than 
the duration of the interstellar shocks (e.g., galactic spi¬ 
ral shocks and supernova remnants). The basic properties 
of the thermal instability are studied in a pioneering pa¬ 
per by Field (1965). Schwarz, McCray, and Stein (1972) 
studied numerically the growth of condensation in cool¬ 
ing region including the effects of ionization and recom¬ 
bination. Hennebelle & Perault (1999) studied the ele¬ 
mentary condensation process in turbulent flow in the re¬ 
stricted conditions of neutral atomic gas in plane-parallel 
geometry. Burkert & Lin (2000) studied cooling and frag¬ 


mentation of gas using simplified power-law cooling func¬ 
tion. Smith (1980) studied collisions between cold atomic 
clouds, which produce thick layers of shock-heated atomic 
gas in which thin sheets of cold molecular gas form by the 
thermal instability. Koyama & Inutsuka (2000, hereafter 
Paper I) have done one-dimensional hydrodynamic calcu¬ 
lations for the propagation of a strong shock wave into 
warm neutral medium (WNM) and cold neutral medium 
(CNM) including detailed thermal and chemical processes. 
They have shown that the post-shock region collapses into 
a cold layer as a result of the thermal instability. They ex¬ 
pect that this layer will break up into very small cloudlets 
which have different translational velocities. 

In this paper, we show that the fragmentation of the 
shock-compressed layer indeed provides turbulent conden¬ 
sations, by using two-dimensional hydrodynamic calcula¬ 
tion with radiative cooling and heating and thermal con¬ 
duction. 


2. NUMERICAL SIMULATIONS 


2 . 1 . Numerical Scheme 


The hydrodynamics module of our scheme is based on 
the second-order Godunov method (van Leer 1979). We 
solve following hydrodynamic equations. 


g + V>v) = 0, 

<9v VP 

— + v • Vv =-, 

at p 


(1) 

( 2 ) 


^+V-(pev)+PV-v = -£-r- ( A(T)+v-(iwr), 

at mu V m H/ 

(3) 

where p, P, v are the density, pressure, and velocity of 
the gas, the specific internal energy e = P/(y — l)p, 
with 7 = 5/3. mu is the hydrogen mass in gram. T 
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and A are heating and cooling rate, respectively. For 
the coefficient of thermal conductivity we adopt K = 
2.5 x 10 3 T 1 / 2 erg cm - 1 K - 1 s -1 (Parker 1953). Effects of 
self-gravity and magnetic field are not treated in this pa¬ 
per. 

In Paper I, we included the following processes in the hy¬ 
drodynamic calculations: photoelectric heating from small 
grains and PAHs, heating and ionization by cosmic rays 
and X-rays, heating by H 2 formation and destruction, 
atomic line cooling from hydrogen Lycq C II, O I, Fe II, 
and Si II, rovibrational line cooling from H 2 and CO, and 
atomic and molecular collisions with grains. However the 
following analytic fitting function adopted here reproduces 
the features relevant for our purpose: 

A(T)/r = 10 7 exp (-114800/(T+ 1000)) 

+ 14v / Texp (—92/T) cm 3 , (4) 

T = 2 x 10 - 26 ergs/s, (5) 

where T is in Kelvins. 

We use 2048 x 512 Cartesian grid points covering a 1.44 
x 0.36 pc region so that the spatial resolution is 0.0007 pc 
= 140 AU. 

We have tested the linear growth of the thermal insta¬ 
bility (Field 1965). The test is performed with an eigen¬ 
function by using 16 grid points. The numerical results 
reproduce the linear growth rate with relative errors of 
0.05 % when Sp/p < 0.01 

2.2. Initial and Boundary Conditions 

To investigate the shock propagation into WNM we con¬ 
sider a plane-parallel shocked layer of gas. We take the 
y- axis to be perpendicular to the layer. An uniform flow 
approaches the layer from the upper side with a velocity 
V y = —26 km/s and density n = 0.6 cm -3 and tempera¬ 
ture T = 6000 K. We assume that the lower side of the 
shocked layer is occupied by hot, tenuous gas (density n = 
0.12 cm -3 ) with a high pressure P/k b = 4x 10 4 Kcm -3 . 

We set up density fluctuations as follows: 

^Alll = -f— jr [p + kf) T sin(a:fci + ykj + 0(ki, fcj)), 

po %ax i j—o 

(6) 

where k\ = 2m/L is a wave number, 0 is a random phase, 
n is a spectral index, and L is a calculation box size of x- 
direction. We imposed out-going boundary condition for 
the y- direction, and periodic boundary condition for the 
x- direction. 

2.3. Results of Thermal Instability and Fragmentation 

Figure 1 shows density distribution at t =0.337, 1.06 
Myr. We used an initial fluctuation A = 0.05, the power 
index n = 0, and i max = 16. Thermally collapsed layer 
and its fragments are seen in the shock-compressed layers. 

Figure la shows the beginning of the formation of cold 
thermally collapsed layer. Density discontinuity at y ~ 
0.46 pc corresponds to the shock front. Density disconti¬ 
nuity at y ~ — 0.1 pc is the contact surface between the 
shock-compressed layer and the hot interior. Behind the 
shock front, cooling dominates heating and temperature 
decreases monotonically. The ISM in the range of 300 - 


6000 K is thermally unstable. Thus, thermally collaps¬ 
ing layer is subject to thermal instability. The conden¬ 
sations continue to collapse and cool until the radiative 
cooling balances the radiative heating, in such a way that 
approximate pressure equilibrium is maintained between 
condensations and the surrounding gas. Temperature of 
these clumps finally reaches the value of thermal equilib¬ 
rium (T = 20 K at n = 2000cm -3 ). When the cooling 
balances the heating, the condensation ceases to collapse 
but the surrounding warm gas continues to cool radiatively 
and accrete onto the cloud cores. 

The cold clumps have considerable translational veloc¬ 
ity dispersion. The typical velocity dispersion is about 
several km/s. The linear analysis predicts that perturba¬ 
tion with sufficiently long wavelength grows exponentially. 
However, the velocity of the non-linearly developed per¬ 
turbation has the upper limit that is the sound speed of 
the warmer medium (« lOkm/s), because the driving force 
of the instability is the pressure of the less dense warmer 
medium. This non-linearly developed perturbation pro¬ 
duces the translational velocities of the cold clumps. These 
velocity on the order of the sound speed of the WNM is 
highly supersonic with respect to the sound speed of cold 
gas. Thus we can understand why the supersonic velocity 
dispersion of the cold medium is comparable to the sound 
speed of the WNM. 

The cold cloud cools the interface via thermal conduc¬ 
tion. The pressure of the warm medium between two cold 
clumps becomes smaller, which in effect provides an at¬ 
tractive force of two adjacent cold clumps. This attractive 
force is the driving force of the coalescence of two cold 
clumps. Thus, cold clumps coalesce into larger clouds. 

The sizes of these clumps are a few orders of 
magnitude smaller than the Jeans length, Aj = 
1.2 pc^(T/20K)(2000cm - 3 /n) . Thus, each cloud is grav¬ 
itationally stable. The shock front itself remains stable at 
least in 10 6 years in this plane-parallel problem. 

For comparison, we changed the initial fluctuations in 
the range of the amplitude (A=0.05-0.005), and the power 
index (n= 0 ,=b 2 ), but the resultant velocity dispersion de¬ 
fined by f pv^mxdxdy/ f pdxdy changed only by 10 %. We 
also changed the spatial resolution (512x128, 1024x256, 
and 2048x512) but the velocity dispersion changed only 
by 10 %. Thus, Figure 1 shows a typical result of this 
problem. From these results we can demonstrate that the 
role of thermal instability behind the shock front is one of 
the most important processes in the evolution of interstel¬ 
lar shocked layers. 

3. DISCUSSION 
3.1. Formation of Molecules 

The main coolant in cold dense clouds is supposed to be 
CO (see, e.g., Neufeld, Lepp, & Melnick 1995; Paper I). 
H 2 enables the formation of CO and other molecules. The 
abundance of H 2 in the CNM critically depends on the op¬ 
tical depth and local density. The typical column density 
of the layer becomes 2 x 10 2 °cm -2 . We calculated detailed 
thermal-chemical equilibrium state in the dense region by 
adopting the values of density, temperature, and column 
density from this dynamical simulation. Twenty-five per¬ 
cent of the hydrogen is H 2 and 0.03% of the carbon is CO 
in thermal and chemical equilibrium. 
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We simulate the observation of molecular clumps at 
t = 1.06 Myr (Figure lb) as the Position-Velocity (P-V) 
diagram of the 12 CO J = 1 - 0 emission (Figure 2a). We 
use the usual relationship Tp = |A 2 //&b, with the specific 
intensity / = n 2 AL/(47rAv), where n 2 A is the cooling per 
unit volume, L is the path length, Av is the Doppler line 
width. We adopt Av to be the sound velocity. As shown in 
Figure 2a, the typical velocity dispersion in the clumps is 
about a few km/s. The cloud to cloud velocity dispersion 
is also a few km/s. Figure 2b shows that the P-V diagram 
of hydrogen nuclei. The velocity dispersion of the diffuse 
warm medium is about 10 km/s. 

3.2. An Origin of “Turbulence” in Interstellar Clouds 

Diffuse ISM in the Galaxy is frequently compressed by 
supernova explosions (McKee & Ostriker 1977). There¬ 
fore, the shock propagation into the ISM plays an im¬ 
portant role in the evolution of the Galactic ISM. We 
have studied the shock propagation into WNM by two- 
dimensional hydrodynamic calculations, and have shown 
that the thermally collapsed layer breaks up into very 
small cloudlets. Fragmentation of the thermally collaps¬ 
ing layer is a result of the thermal instability. The ther¬ 
mal instability produces many cloudlets which have differ¬ 
ent translational velocities. We expect that the Galactic 
ISM is occupied by these small cloudlets which have super¬ 
sonic velocity dispersion , because the ISM is frequently 
compressed by supernova explosions, stellar winds, spiral 
density waves, cloud-cloud collisions, etc. The H I 21-cm 
observation maps show the existence of many shell-like or 
filamentary structures in the Galaxy (Hartmann & Burton 
1997). In addition, high resolution observations suggested 
that the clumpy distribution is ubiquitous in the Galaxy 
(Heiles 1997). If these structures are really the results of 
the shock waves, many small molecular cloudlets should 
be hidden in the shock-compressed layers. These cloudlets 
should have translational velocities owing to the fragmen¬ 
tation of thermal instability. Observational “turbulence” 
in the ISM should reflect these motions. Thus, we propose 


that an origin of “turbulence” in the ISM is the motion of 
this small cloudlet complex. 

In the radiative shocked layers, the gases lose thermal 
energy through radiative cooling. Thus, the initial kinetic 
energy of the pre-shock gas (in the comoving frame of the 
post-shock gas) converted to the radiation energy, which 
escapes from the system. If, however, the post-shock gas 
becomes dynamically unstable by the thermal instability 
as in this paper, some portion of the thermal energy is 
transformed into the kinetic energy of the translational 
motions of the cold cloudlets, which does not easily escape 
from the system. Therefore we can consider the origin 
of interstellar turbulence is due to the conversion of the 
gas energy in supersonic motion. In principle, the kinetic 
energy of the cloudlets can be lost via coalescence of the 
cloudlets. The damping rate of the velocity dispersion of 
the cloudlets will play a key role in the evolution of this 
system, and hence, should be studied in the subsequent 
paper. 

Expected effects of magnetic field modify the develop¬ 
ment of the thermal instability described here. The pres¬ 
ence of magnetic field can suppress thermal conduction 
efficiently, allowing the collapse of small scale structure by 
the thermal instability. The magnetic forces suppress the 
motion that is perpendicular to the magnetic field lines, 
and hence suppress the growth of perturbations whose 
wavevector is perpendicular to the magnetic field lines. 
However, the magnetic forces do not affect the motion that 
is parallel to the magnetic field lines, and hence the growth 
of perturbations whose wavevector is parallel to the mag¬ 
netic field lines. As a result, the dense sheets or filaments 
will form and tend to align perpendicular to the magnetic 
field lines. We need three-dimensional calculation to ana¬ 
lyze these effects of magnetic fields. 

Numerical computations were carried out on 
VPP300/16R and VPP5000/48 at the Astronomical Data 
Analysis Center of the National Astronomical Observa¬ 
tory, Japan. 
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Fig. 1. — Density distribution of the shock-compressed layers. The figures show the snapshot at (a) t = 0.337 Myr, (b) t = 1.06 Myr, 
respectively. 



Fig. 2.— The panel a shows the position-velocity (P-V) diagram of the 12 CO J=1 - 0 emission obtained from the calculation at t = 1.06 
Myr (see Figure lb ). We assume 0.03 % of the carbon is CO. The panel b shows the P-V diagram of the hydrogen nuclei. 


























